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i-rt I We numerically investigate the problem of the propagation of a shock 

^ • in an horizontal non-loaded granular chain with a bead interaction force 

exponent varying from unity to large values. When a is close to unity 
we observed a cross-over between a nonlinearity-dominated regime and a 
solitonic one, the latest being the final steady state of the propagating 
^ . wave. In the case of large values of a the deformation field given by the 

£^ I numerical simulations is completely different from the one obtained by an- 

alytical calculation. In the following we studied the interaction of these 
shock waves with a mass impurity placed in the bead chain. Two different 
, physical pictures emerge whether we consider a light or a heavy impurity 

Q>^ I mass. The scatter of the shock wave with a light impurity yields damped 

oscillations of the impurity which then behave as a solitary wave source. 
Differently an heavy impurity is just shifted by the shock and the trans- 
mitted wave loses its solitonic character being fragmented into waves of 

'^ ■ decreasing amplitudes. 

o 
o 



1 Introduction 

rS \ Since the pioneering work of Nesterenko |jl| there has been an increase of activity 

c^ ' devoted to the study of wave propagation in granular media. The theoretical 

prediction of solitons-like waves propagating in one dimensional chains of perfect 
elastic spheres has been experimentally checked by Lazaridi and Nesterenko [^ 
and more recently by Coste, Falcon and Fauve 0. These studies concern the 
idealized case of spheres interacting according to the Hertz law where the inter- 
action force depends on the overlap distance of the spheres to the power |. The 
result of Nesterenko for the particular Hertzian case can actually be generalized 
to an interaction law where the exponent can vary from one to infinity [0], the 



case of the exponent equal to unity being special since it leads to a self-similar 
asymmetric pulse propagation resulting from a shock perturbation |^, ||. Instead 
of traveling as a soliton the pulse propagates increasing its width and decreasing 
its amplitude. It was then showed that the amplitude decreases like t~°'^^ and 
the width increases like ts. This generalization of the interaction law by letting 
the interaction exponent vary allows us to take into account the various shapes 
of two interacting beads. 

Besides isolating non-linearities in one dimensional chains other physical ingre- 
dients have been added in order to make these models more realistic. Sinkovits 
and Sen |^ did molecular dynamics simulations on gravitationally compacted 
Hertzian chains. It was then observed that gravitation makes waves loose their 
soliton-like features becoming solitary waves with the wave speed increasing with 
depth z as ze as predicted by Hertz law. Besides, if an impurity is placed in the 
chain, the propagating wave is scattered by it and a reflected wave is obtained 
in addition to the transmitted wave, these waves having a similar shape as the 
incident wave. Following these first results Hong, Ji and Kim |^ extracted several 
power laws for the velocity amplitude of a wave propagating in a gravitationally 
compacted chain. It was shown that in a homogeneous chain of beads the ampli- 
tude velocity of the solitary pulse decreases like 2:~4(3 + a) for an interaction law 
exponent a greater than unity. This power law is only valid in the case of a small 
perturbation where the non-linear interaction law can be locally linearized, the 
stiffness coefficient then varying with depth. 

In this paper, we present molecular dynamics results for non-loaded horizontal 
granular chains containing impurities. In our model the force interaction law ex- 
ponent is allowed to vary from one to infinity, the Hertzian case being considered 
as a special one. After having reviewed the theoretical results for shock propaga- 
tion in a horizontal bead chain and their relevance in reproducing the results from 
numerical simulations we study the interaction of shock waves with impurities of 
various masses. We show that the patterns observed for the scattered waves vary 
as impurity masses are larger or smaller than the mass of the beads constituting 
the granular chain. In the case of a light impurity we observed a backscattered 
solitary wave followed by an increasing number of secondary waves of decreas- 
ing amplitudes produced at the impurity location. These secondary waves also 
propagate in the forward direction following the transmitted wave. We looked 
at the displacement of the impurity after the collision. The impurity actually 
oscillates in a ballistic motion around its equilibrium position. These oscillations 
are damped since the impurity collides its nearest neighbors producing a back- 
ward and a forward moving pulse in each cycle. If we consider the interaction of 
a shock with an heavy impurity we obtain a very different scenario. Instead of 
exciting the impurity the shock only shifts it in its moving direction. The chain 
is then halved in two separated chains on which propagate several pulses. The 
situation is no more symmetric as in the previous case. After the collision the 
incident wave does not produce a stable transmitted wave but is decomposed in 



a series of pulses of decreasing amplitudes. Before the impurity we observe the 
propagation of a stable reflected solitary wave. In each case the amplitude of the 
reflected pulse grows with the absolute value of the difference between the mass 
of the impurity and the normal mass. Hence the speed of the reflected pulse 
grows with this mass difference since the speed of our solitary waves increases 
with the amphtude. 

In section 2 we present the model and its analytical solution in the continuum 
limit. This solution is compared to the numerical simulations as the interaction 
law exponent is tuned from unity to large values. In section 3 we analyze the 
interaction of these shock waves with an impurity and display the energy distri- 
bution in the bead chain after the scattering with the impurity. A conclusion is 
given in section 4. 

2 Theoretical results versus numerical simula- 
tions 

The horizontal granular chain consists of N touching beads of equal sizes and 
masses. The repulsive interaction force between two adjacent beads under com- 
pression is written F = kS"" where 6 is the overlap of the beads. For the exponent 
a equal to | we have the Hertz law for perfectly spherical beads and the stiff- 
ness constant k can be expressed as a function of the spheres radii, the Young 
modulus and the Poisson ratio of the material. Here we let a vary from one to 
infinity considering beads with different geometries. The equation of motion for 
the displacement Ui of the bead i is given by 

miii = k5tiA^:-i,{) - fc<^,Vi^('^M+i) < z < TV - 1 (1) 

In terms of the displacements of the beads i — \ and z + 1 one can write 

mui = k{ui^i - Ui)"9{ui^i - Ui) - k{ui - Ui+i)°'9{ui - Ui+i) (2) 

the Heaviside function having been introduced in order to have no tension forces. 
This equation has been solved by Nesterenko [^ in the case of a shock per- 
turbation by doing a long wavelength development. This approximation will be 
valid if the spatial extension of the wave generated by the shock is much greater 
than the size of a single bead. Assuming this we write Ui as u{x,t). Moreover, 
since we are solving the equation for a wave in which the beads are compressed 
one can get rid of the Heaviside functions. The resulting equation is written 
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where a is the dimension of the bead and we developed u{x + a) and u{x — a) to 
fourth order. We then rewrite the equation 
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Ux being of order ^ where L is the dimension of the pulse and expand the powers 
in Taylor series to the order ■^. We then obtain 
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Looking for traveling wave solutions u{^) = u{x — Vt) with V being the speed of 
the wave, a longer calculation [Q gives us for the deformation 
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This periodic wave solution is singular since it violates the condition —Ux > 
in eq.(P). We actually get the shock wave solution by imposing the argument 
of the cosines to vary in the interval [— |(2fc + 1), |(2A; + 1)], k E Z. From the 
expression above we easily extract the wave speed V: 
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Vmax being the maximum of the velocity pulse. 

The accuracy of this analytical approximation has been tested numerically 
with a Gear fifth order predictor-corrector algorithm on a horizontal chain of 200 
beads with free boundary conditions. A shock perturbation is imposed on the 
chain by giving the first bead a 5m. s'^ speed at zero time. We choose beads of 
mass m = Ig and stiffness k = lO^A^.m"^. We first looked at the shape of the 
deformation profile of the bead number 100 as function of time and compared 



it with with the shape given by the expression (|^ for different values of a. 
We displayed in fig.(|l]), fig-(|D and fig.(@) the superposition of numerical and 
analytical curves for a equal to 1.2, 1.5 and 1.7. One easily sees on these curves 
that the analytical results are rather good for low values of a (a > 1) showing 
discrepancies when a is increased. Actually, the tails of the deformations curves 
are much longer than the ones given by equation (0). If the value of a is increased 
to higher values the numerical and analytical results become completely different 
as one can see in fig.(P where we took a equal to 10.0. Other discrepancies can 
be found if a is taken in the neighborhood of unity. It can be seen in fig.® 
that the deformation curve on the bead number 100 is not symmetric when a is 
taken equal to 1.01. This non-symmetric shape actually evolves to a symmetric 
one as the pulse propagate. We simulated a chain of 3000 beads and considered 
the deformation of the bead number 2000 : the pulse has the soliton shape given 
by equation (|^) (fig.(^). The pulse has actually to propagate over 2000 beads 
before reaching its solitonic regime predicted by (^. In fact, there is a cross-over 
between a nonlinearity-dominated regime and a solitonic regime where dispersion 
is balanced by nonlinearity. The non-linearity dominated regime corresponds to 
a transient where the amplitude of the pulse decreases slowly in time before 
reaching a constant value, i.e, the solitonic regime. At the same time the width 
of the pulse increases also reaching a stationary value. Moreover the shape of the 
pulse evolves continuously from a non-symmetric one to a symmetric one. The 
time of the transient obviously increases as we decrease a (if we take a = 1.004 
the pulse has to propagate over 15000 beads before reaching its stationary value). 
This cross-over behavior can be understood by looking at equation (H) for a close 
to one. We can rewrite this equation with a small parameter e where a = 1 + e : 
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If we consider {—UxY constant equal to one we obtain the long wavelength equa- 
tion for the linear springs chain which is the same as the simple wave equation 
up to a dispersive term, the fourth order derivate of the displacement. It has 
been shown |^, |^ that the velocity solution of this equation is not a soliton but 
a non-symmetric self-similar wave of decreasing amplitude and increasing width, 
the amplitude going to zero and the width going to infinity as time goes to infin- 
ity. For a close to one the non-linearity {—UxY vary slowly and the time for it to 
balance the action of the dispersive four order term increases as e goes to zero. 
We can actually calculate how the amplitude and the width of the velocity pulse 
vary with epsilon. The velocity of a bead is given by eq.(0) up to the factor V . 
Replacing a with e + 1 we get : 
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where 

-(^£)^ 

We then write eq. ([TID in an exponential form and expand its argument in a Taylor 
series keeping the lower order term : 

The solution of eq.(||) is then approximated by a Gaussian for a close to unity. 
According to the properties of Gaussian distributions the width w of the solitary 
wave follows the scaling law w oc e~2, i.e : 

w(x{a-l)~^ (14) 

For a = 1 we recover the result of [^ Q where the width of the wave was found 
to diverge as time goes to infinity. We checked numerically the accuracy of this 
scahng relation and found a good agreement (see Fig.(|^). A second scaling 
relation can be obtained by using the conservation of the kinetic energy in the 
wave since the beads following the wave have negligible velocities. It follows that 
the product v'^^^w remains constant. Hence, 

Vmax OC {a -1)^ (15) 

We again recover the result of P, § where the velocity amplitude was found to 
go to zero as time goes to infinity. We displayed in Fig. (|^) the value of the final 
velocity amplitude Vmax when a is varied in the neighborhood of one and found 
again a good agreement with the obtained scaling relation. These two scaling laws 
show that the transition between the self-similar regime a = 1 to the solitonic 
regime a > 1 is actually continuous, the value a = 1 being not singular. 

3 Collision of a shock wave with an impurity 

After having considered the propagation of a single solitary wave in a chain of 
equal masses we want to present what are the consequences of the introduction of 
a single impurity on the stability of a shock wave. For that, we consider a chain 
of 1000 beads, the bead 500 being replaced by an impurity. By impurity we mean 
a bead whose mass is larger or smaller than the mass of its neighbors. Different 
behavior arises as we introduce a light or an heavy impurity. Besides, we also have 
to take into account the value a of the interaction law exponent. The simplest 
way of looking at what happens when a shock generated on the left end of the 
chain scatters the impurity is to compute the velocity and the displacement of the 



beads after the collision. Fig.(^ displays velocities for a light impurity and an 
exponent a equal to 1.5. The leading wave corresponds to the shock wave which 
is followed by an increasing number of secondary solitary waves with decreasing 
amplitudes generated at the impurity position. A symmetric pattern arises in 
the opposite direction since solitary waves with negative velocities and almost 
same amplitudes as the previous ones propagate. The wave speed of the incident 
shock is modified after the collision since its amplitude decreases providing with 
energy the impurity which can then generate new waves. The corresponding 
displacement profile is displayed on Fig.(0) where the pulses are now replaced 
by the squared front waves characteristics of shock waves. In order to better 
understand the phenomenon we looked at the spatio-temporal pattern exhibited 
by the scattering. In Fig. ( pH]) one can see black straight lines defining compression 
waves whereas white zones mean beads are at touching. The secondary pulses are 
created at the impurity position and after a transient displayed by thin curved 
black lines near the impurity location they propagate with a constant velocity 
reaching then a stationary regime. The wave speed of these pulses are decreasing 
as new pulses are created as one can see by looking at the slopes of the black lines. 
It agrees with the display of waves with decreasing amplitudes in Fig.(^. From 
these observations we can conclude that the light impurity acts as a wave source. 
In order to understand the behavior of this source we computed the displacement 
of the impurity as function of time after the scattering. Fig.(|12|) clearly shows 
that the impurity oscillates around an equilibrium position. These oscillations 
are zig-zag shaped showing that the impurity moves with a constant velocity 
between two shocks with its left and right neighbors. This velocity decreases in 
time (lower slopes of the oscillations) since the impurity transfers energy to its 
neighbors at each shock creating new pulses. Although we chose special values 
for the impurity mass and the interaction law exponent the same behavior is 
observed if we maintain the impurity mass lower than the chain masses and the 
interaction exponent greater than one. 

If we now turn to the case of an heavy impurity we observe a very different 



physical picture. In Fig.(|T3D we displayed the velocity of the beads for an impurity 
mass greater than beads chain mass, the interaction exponent remaining equal 
to 1.5. The collision generates many pulses of decreasing amplitudes behind the 
impurity, these pulses being very close to each other. Moving in the opposite 
direction one can see one pulse of large amplitude. This pattern is not symmetric 
as in the case of a light bead and a different physical process seems to be involved. 
We can better understand what happened by looking at the spatio-temporal 
pattern of contacts displayed in Fig.fpJlf). After collision the bead chain is halved 
as it is shown by the white vertical line separating two black areas meaning 
compression between beads. The impurity has actually been given a shift to the 
right when it collided with the shock wave. Instead of bouncing back the impurity 
sticks and there is the formation of a compression chain. The beads following the 
impurity don't remain homogeneously compressed and there is a fragmentation 



of the incident pulse giving raise to many pulses corresponding to the straight 
lines of decreasing slopes on the spatio-temporal pattern. These lines obviously 
correspond to pulses of decreasing amplitude in Fig. ([T3|) . Between these pulses 
the beads remain at rest just touching their nearest neighbor. Other compression 
areas appear on the left side. The straight line of larger slope is the main reflected 
wave. A growing compression zone moving at constant velocity follows this pulse. 
After some times this area also begins to fragment into new pulses as one can 
see from the black straight lines appearing at the front. These new compression 
pulses are not visible on the velocity profile since their amplitude is too low. This 
spatio-temporal pattern allows us to say that the collision of a shock wave with an 
heavy impurity make it fragment into solitary waves of decreasing amplitudes. 
This is not what happens when we introduce a light impurity since the shock 
excites the impurity. To our knowledge, the scattering with an heavy mass is 
the simplest way to make the solitary wave unstable. As in the previous case, 
although we took lOg for the impurity mass and 1.5 for the interaction exponent 
the obtained pattern remains unchanged as these two parameters are tuned. 
The number of fragmented pulses actually increases with the mass of the impurity. 
To show it we considered a chain of equal masses with a bead of larger mass at 
its beginning. In the case of a bead having the same mass we know that the 
perturbation of the first bead by giving it a finite velocity at zero time will 
produce a single pulse. As we increase the mass of this first bead an increasing 
number of pulses appear by the same process as described previously. We plotted 
the corresponding velocity profiles for two different values of the first bead mass 
in Fig.(|l|) and Fig.®. 

We also looked at the energy distribution in the bead chain as function of the 
mass of the impurity and of the value of the force interaction exponent. Fig.(0) 
shows the ratio of the reflected energy at the impurity location over the total 
energy in the chain. The reflected energy is defined as the energy contained in 
the first half of the bead chain after scattering. This energy is not stationary 
immediately after wave reflection and we have to wait the time for it to reach 
a stationary value. The impurity mass varies between 0.6g and 5g and a takes 
values between 1.01 and 5. One easily sees that the amount of reflected energy 
increases as the mass difference increases. This amount also increases with a if 
the the impurity mass is greater than a value located between 1.5g and 2g. Below 
this value the curves cross each others at different values of the impurity mass. 
Below the crossing it is the lowest value of the two exponents which gives the 
greater amount of reflected energy. If we now consider the ratio of reflected energy 



over transmitted energy in Fig. ([T8|) we observe similar variations as previously. 
Moreover the ratios remain lower than unity meaning that the amount of reflected 
energy is always lower than the amount of transmitted energy in the impurity 
mass interval considered. 



4 Conclusion 

We investigated the problem of a shock propagating in a non-loaded horizontal 
granular chain with a general beads interaction force F oc 5", a > 1. We com- 
pared our numerical simulations with previous analytical results obtained in the 
long wavelength approximation [Q. For values of a close to unity we found a 
cross-over between a non-linearity dominated regime and a solitonic one corre- 
sponding to the long wavelength approximation. The latter regime is reached 
after a time which increases and tends to infinity as one tunes a to unity. It was 
then shown that the transition between the a = 1 P, |^ (self-similar) regime to 
the a > 1 regime is continuous. In the case of large values of a we observed dis- 
crepancies between numerical simulations and analytical results. As a increases 
beyond 1.2 tails in the deformation curves separate. When a reaches the value of 
1.7 the analytical result given for the deformation curve is already not reliable. 
In the limit of very large values of a the deformation curve takes a triangular 
shape totally inconsistent with the one given by the theoretical analysis. 
We then studied the interaction of these solitonic shock waves with an impurity 
mass placed in the middle of the bead chain. We first considered light impurities. 
The scattering of the shock with the impurity actually excites it in an oscillatory 
ballistic damped motion. The impurity acts as a wave source emitting solitary 
waves of decreasing amplitude in both directions. If we place an heavy impu- 
rity in the bead chain we observe a completely different pattern. The impurity 
is no more excited but sticks its next neighbor splitting the bead chain in two 
parts. The chain corresponding to the transmitted wave supports no more soli- 
tonic propagation and we have an instability displayed by the fragmentation of 
the transmitted wave into waves of decreasing amplitudes. On the other half of 
the chain a reflected stable solitonic wave propagates followed by a growing un- 
stable fragmented front appearing at the impurity location. The physics involved 
in this case is much more complicated and deserves further investigations. 
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Figure 1: Deformation curves of the bead number 100 in a chain of 200 beads, 
the force exponent being equal to 1.2. There is very good agreement between 
numerical and theoretical results. The crossed curve from the analytical result 
eq.(|^) superimposes the plain line curve from the numerical simulation. The 
crossed curve obviously correspond to a single oscillation of the periodic solution 
eq.(0). 
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Figure 2: Deformation curves when a = 1.5. Discrepancies between numeri- 
cal and theoretical results appear at the curves tails. The crossed tails of the 
analytical result remain below the real tails of the numerical simulation. 
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Figure 3: Deformation curves when a = 1.7. The increase in a increases the 
deviation of the curves at the tails. The oscillating character of the analytical 
result is displayed at the corners of the plot where one can see the end of the 
previous oscillation of the non-linear periodic solution on left and the beginning 
of its next oscillation on right. 
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Figure 4: Deformation curves when a = 10. There is no agreement between 
analytical and numerical results. The true deformation has a triangular shape 
while the analytical curve display a periodic parabola shape. 
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Figure 5: Deformation curve when a = 1.01 for bead number 100. The pulse is 
not symmetric. We are in the non-hnearity dominated regime. 
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Figure 6: Deformation curves when a = 1.01 for bead number 2000 of a bead 
chain of 3000 beads. The pulse reached its stationary sohton shape given by 
eq.(0) and the numerical and analytical curves superimpose. 
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Figure 7: Plot of the width of the pulse as function of a when a is close to unity. 
The dashed line displays the corresponding power law : w (x[a — i) -0.492^ 
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Figure 8: Plot of the velocity pulse amplitude as function of a when a is close to 
unity. The dashed line displays the corresponding power law : v^ax oc (a — 1)°'^^^. 
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Figure 9: Velocities of the 1000 beads after the scattering of the shock with a 
light impurity of mass rrii = 0.6g, the normal mass being m = Ig. The force 
interaction exponent is a = 1.5. Sohtons are emitted at the impurity location in 
forward and backward directions. 
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Figure 10: Displacement profile corresponding to the velocity profile. 
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Figure 11: Spatio-temporal evolution of beads contacts for the collision with a 
light impurity of mass rrii = 0.6g. A black dot means that two beads are in 
compression while a white dot means that they are at touching, a is equal to 
1.5. Black straight lines correspond to compression pulses. 
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Figure 12: Displacement of the impurity located at the position 500 in function 
of time for rrii = 0.6g and a = 1.5. The light impurity oscillates in a damped 
ballistic motion showed by the zig-zag shape of the curve. 
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Figure 13: Velocities of the 1000 beads after the scattering of the shock with 
a heavy impurity of mass rrii = lOg {a = 1.5). The transmitted sohtary wave 
fragments into many sohtary waves of decreasing amphtudes. The major reflected 
part is displayed by a single pulse. 
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Figure 14: Spatio-temporal evolution of beads contacts for the collision with an 
heavy impurity of mass nii = lOg {a = 1.5). In addition to the transmitted 
fragmented front on can see a second fragmenting front following the reflected 
wave. 
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Figure 15: Velocities in a 100 beads chain when it is perturbed with an heavy mass 
of 3g {a = 1.5) with initial velocity of 5r?7,.s~^. Instead of having a single soliton 
as it is the case for a homogeneous chain many pulses of decreasing amplitudes 
are propagating. 



25 




10 20 30 40 50 60 70 80 90 100 



Figure 16: Velocities in a 100 beads chain when it is perturbed with an heavy 
mass of lOg {a = 1.5) with initial velocity of 5m. s^^. As the mass of the first 
bead is increased the number of emitted pulses grows. 
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Figure 17: Ratios of the reflected energy over tlie total energy of the bead chain 
for different values of a. Going from the upper right curve to the lower right 
curve a takes the values 5, 3, 2, 1.5, 1.2, 1.1 and 1.01. When mj is greater than a 
value between 1.5 and 2 the energy ratio increases with a. Below this threshold 
we have a more complicated pattern in which curves cross each others at different 
locations. 
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Figure 18: Ratios of the reflected energy over the transmitted energy. Going from 
the upper right curve to the lower right curve a takes the values 5, 3, 2, 1.5, 1.2, 
1.1 and 1.01. A similar pattern as the one for the reflected energy over the total 
energy is obtained. The part of reflected energy remains lower than the part of 
transmitted energy in the impurity mass considered. 
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